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The interplay between topology and dynamics in complex networks is a fundamental but poorly 
explored problem. Here we study this phenomenon on a coupled model describing extremal dynamics 
on a network which is in turn shaped by the dynamical variable itself. Each vertex is assigned a 
fitness, and the vertex with minimum fitness and its neighbours are updated as in the Bak-Sneppen 
model. Simultaneously, the fitness values determine the connection probability as in the fitness 
network model. We solve the model analytically and show that the system self-organizes to a 
nontrivial state which differs from what is obtained when the two processes are decoupled. A 
power-law decay of dynamical and topological quantities above a threshold emerges spontaneously, 
as well as a feedback between different dynamical regimes and the underlying network's correlation 
and percolation properties. Our model is a prototype for more general mechanisms exploring such 
unexpected effects in networks. 



The properties of dynamical processes defined on com- 
plex networks display a strong dependence on the topol- 
ogy [E S H, I3l- On the other hand, there is growing 
empirical evidence S 0] that many networks are in 
turn shaped by some variable associated to each vertex, 
an aspect captured by the 'fitness' or 'hidden- variable' 
model [1, Q. Until now, these two facets of the same 
problem have been treated as separate, by considering on 
one hand dynamical processes on static networks [l|, |j| , 
and on the other hand network formation mechanisms 
driven by quenched variables d, d, [13, [HI [13] ■ This may 
be perhaps justified for short time scales. However, in 
the long-term evolution it is crucial to understand the 
effects that these mechanisms have on each other, with- 
out ad hoc specifications of any fixed structure either in 
the topology or in the dynamical variables. Remarkably, 
the interplay of dynamics and topology can drive the net- 
work to a self-organized state that cannot be inferred by 
studying the two evolutionary processes as decoupled. 

Here we explore explicitly the possibility that the net- 
work supports a dynamical process which in turn shapes 
its topology, with a continuous feedback between dy- 
namics and structure. Models where both dynamical 
and topological properties are continuously updated have 
been considered ,^13, 14, 15, 16, J/7„dS.]. In these cases, 
however, the rewiring of links is not completely driven by 
the dynamical variables. By contrast, our main interest 
here is the description of a self-organized process where 
the dynamical variable fully acts also as the 'hidden vari- 
able' shaping network topology explicitly, as in the fitness 
model. Due to the increased complexity of the problem, 
in this paper we choose the simplest possible dynamical 
rule for the hidden variable. We focus on the extremal 
dynamics defined in the Bak-Sneppen (BS) model llEji. a 
traditional model of self-organized criticality (SOC) [20| 
inspired by biological evolution. Since the outcomes of 
this model on a wide r ang e of fitness-independent net- 
works are well known [11 [HI, [H, [H, [H, [ll, [l^], it is 



straightforward to understand what are the novel effects 
originating uniquely by the interplay with the fitness- 
driven topological evolution we consider here. 



Coupling the Bak— Sneppen model and the Fitness 
network model 

In the traditional Bak-Sneppen model [l3] defined on 
a generic graph [HJ, [H, [H, [3, [H, [11], each of the N 
vertices is regarded as a biological species having a fit- 
ness value Xi , initially drawn from a uniform distribution 
between and 1. At each timestep the species with low- 
est fitness and all its neighbours undergo a mutation, 
and their fitness values are drawn anew from the same 
uniform distribution. The process is iterated, and even- 
tually the system reaches a stationary state, character- 
ized by a step-like fitness distribution, uniform above 
a threshold r. This behaviour is observed on regular 
lattices [H, l2l| r andom graphs [2^ . small-world [13] 
and scale-free [2^ [25l . [26| networks, the only depen- 
dence on the particular topology being in the value of 
T[ii, [m, [13, [ii, [13, [H, [Ig. In particular, r vanishes 
for scale-free degree distributions with diverging second 
moment [13,[13, .26]. 

Here we couple this dynamical rule with the fitness 
model assumption [3] that the network is formed by draw- 
ing a link between any two vertices i and j with fitness- 
dependent probability f{xi,Xj), thus introducing an in- 
trinsic feedback between dynamics and topology. In this 
way, whenever the fitness Xi of a species i is updated to 
x'^, the links from i to all the other vertices j are drawn 
anew with probability f{x[,Xj). Besides the updates de- 
scribed above, one could also define additional and arbi- 
trary link updating events. In other words, in addition 
to the 'natural' link update occurring between a mutat- 
ing species and all other species, other link updates may 
happen between any two vertices i and j at generic and 
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arbitrarily distributed timesteps. When this occurs, the 
hnk between i and j is updated and drawn anew with 
probabihty f{xi, Xj) where Xi and Xj are the current fit- 
ness vahies of i and j, even if the latter are not involved in 
a mutation event. Remarkably, it is possible to show (see 
Supplementary Information) that the introduction of link 
updating events leaves the system in the same stationary 
state as if they were absent. Therefore our model is very 
general in this respect, and allows for rearrangements of 
ecological interactions on shorter timescales than those 
generated by mutations. In particular, the stationary 
state is the same if the whole network is updated at each 
timestep. In this case storing the information on the 
adjacency matrix among species is unnecessary, and we 
shall exploit this property to achieve fast and very large 
numerical simulations of the model. 

As we show below, the coupling between structure 
and dynamics leads to unexpected results that cannot 
be traced back to any of the two processes taken as sepa- 
rate. Moreover, another important advantage is that the 
main limitations of the two models disappear when they 
are coupled together. A fundamental problem in the BS 
model on static graphs is that, after a mutation, the new 
species always inherits exactly all the links of the previ- 
ous one. This is hard to justify, since it is precisely the 
structure of ecological connections among species which 
is believed to be both the origin and the outcome of 
macroevolution [l^]- Here the fitness-driven link updat- 
ing overcomes this problem. Similarly, the static fitness 
model requires the specification of an ad hoc fitness distri- 
bution that never changes. By contrast, here the fitness 
distribution self-organizes spontaneously to a stationary 
probability density, removing the need of arbitrary speci- 
fications. As we discuss below, a proper interpretation of 
the fitness also allows us to remove the remaining arbi- 
trariness in the choice of f{xi,Xj). However, in order to 
keep our approach as general as possible, we first study 
the model analytically for a generic form of /, and focus 
on particular choices only later. 

The analytical solution of the model for an arbitrary 
linking function f(x,y) can be obtained by focusing on 
the master equation for the fitness distribution p{x) at 
the stationary state (see Supplementary Information). 
We find that the analytical expression for p{x) is 



(1) 




N JJ" f{x, m)dm 



where r is a threshold value determined through the nor- 
malization condition p{x)dx, which reads 



dx 



Jo f{x,m)dm 



(2) 



In the infinite size limit N — > oo, the distribution q{m) 
of the minimum fitness value m = Xmin is uniform be- 
tween and r, while all other values (except possibly a 



vanishing fraction) are above r (see Supplementary In- 
formation). In other words, q{m) — Q(t — m)/T. This 
characterizes the stationary state completely. Once p{x) 
is known, all the expected topological quantities can be 
determined as in the static fitness model [ll|, . For 
instance, the average degree of a vertex with fitness x is 
given by 



k{x)^N / f{x,y)p(v)dy 



(3) 



and the inverse function x{k) can be used to obtain the 
cumulative degree distribution as 



P>(fc) ^ 



fc(i) 



P{k')dk' = p>[xik)] 



(4) 



where p>{x) = p{x')dx' is the cumulative fitness dis- 
tribution. Note that if r is nonzero the fitness distribu- 
tion preserves the discontinuous behaviour displayed on 
static networks 19, 21, 22, 23, 24111, HI]. However, here 
we find the novel feature that p{x) is in general not uni- 
form for X > T. This unexpected result holds for any 
nontrivial choice of f{x,y), and hence for any topology. 
Therefore the effect is not due to the topology itself, but 
to the interplay between the topological evolution and 
the dynamical process entangled with it. Remarkably, 
this feedback alone determines the self-organization of 
the system from a random structure to a complex net- 
work with nontrivial dynamical and topological proper- 
ties. 



Fitness— independent random graphs 

The above analytical solution holds for any form of 
f{x,y). Now we consider possible choices of this func- 
tion. First note that the null choice is f{x,y) — p, the 
network being a random graph. It is nonetheless an in- 
structive simple case, and we briefly discuss it. Moreover, 
this choice is asyrnptotically equivalent to the random 
neighbour variant [2^ of the BS model, the average de- 
gree of each vertex being d = p{N — 1) « pN (we drop 
terms of order 1 /N from now on) . Our analytical results 
read 



p{x) 



{tN)- 
{ptN)- 



X < T 
X > T 



(5) 



and, depending on how p scales with N, eq.Q implies 




(6) 



We note that these three dynamical regimes are tightly 
related to an underlying topological phase transition. As 
p decreases, the whole system splits up into a number of 
smaller subsets or clusters. Such process displays a criti- 
cal behaviour near the threshold Pc ~ 1 /N [J, Q . Below 
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Pc, each node is isolated or linked to a small number of 
peers. Above Pc, a large giant component emerges in- 
cluding a number of nodes of order 0{N), whose fraction 
tends to 1 as p — > 1 . This explains the dynamical regimes 
in eq.®. If pN oo (dense regime), then r — > and 
p{x) is uniform between and 1 as in the initial state, 
since an infinite number {kmin) = pN of fitnesses is con- 
tinuously udpated as on a complete graph. In this case, 
the step-like behaviour is destroyed. If pN = d with 
finite d > 1 (sparse regime), then r remains finite as 
N ^ oo, and this is the case considered in ref. [13] that 
we recover correctly. Finally, if p falls faster than 
the graph is below the percolation threshold (subcriti- 
cal regime): the updates cannot propagate and t — > 1, 
as for N isolated vertices (only the minimum is contin- 
uously updated, which after many timesteps results in 
pushing all fitness values, except the newly replaced one, 
towards 1). Therefore the dynamical transition is rooted 
in an underlying topological phase transition. This pre- 
viously unrecognized property is fundamental and, as we 
show below, is also general. 



state. All other higher-order properties are completely 
random, except for the structural correlations induced by 
the degree sequence [13, • It therefore represents the 
fitness-dependent version of the so-called configuration 
model 

U,m- As we show later on, structural correla- 
tions have an important impact on the dynamics. With 
the above choice, p{x) can be directly computed analyt- 
ically through eq.([T]). However we write it in a differ- 
ent form, which is equivalent when A'^ — > c», in order to 
solve also more complicated integrals involving it later 
on. We use the trick {f{x,m)) ~ f{x,{m)) where the 
angular brackets denote an average over the distribution 
q{m) of the minimum fitness, that is /(x, m)dm « 

{zxt/2)/{1 + zxt/2). As we show in a moment, when 
N ^ oo this approximated expression becomes exact. 
Then eq.((T]) yields 



p{x) 



(rAf)-i X <T 

(tN)-^ + 2/{zNt^x) x>t 



(8) 



where r is the solution of cq.®, which reads 



Fitness— dependent complex networks 

A nontrivial form of f{x,y) must be chosen carefully. 
On static or fitness-independent networks Xi is usually 
interpreted as the fitness barrier against further muta- 
tion, and the links are interpreted as feeding relations 
[l9j . However, once the topology depends on x these two 
interpretations are difficult to concile. The coupling we 
have introduced requires consistent interpretations of x 
and of the links. Also, the form of f{x,y) must be con- 
sistent with the feature that the updates of x propagate 
through the network determined by it. This very instruc- 
tive aspect must characterize any model with coupled 
topology and dynamics, and reduces significantly the ar- 
birariness introduced in the static case. Here we suggest 
that the simplest self-consistent interpretation is the fol- 
lowing. Since there is no external world in the model, the 
environment experienced by a species is simply the set of 
its ecological interactions. Now let Xi represent the fit- 
ness (rather than the barrier) of i, and let a link between 
two species mean 'being fit to coexist with each other' 
(i.e. it represents an undirected, non-feeding interaction 
beneficial to both). The more a species is connected to 
other species, the more it is fit to the environment. This 
picture is self-consistent provided that the larger x and 
y, the larger f{x, y). Following the results of refs.[27l.[28l|. 
the simplest unbiased [111 choice for such a function is 



f{x,y) = 



zxy 
1 -f zxy 



(7) 



where z is a positive parameter controlling the number 
of links. This choice generates a network with a nonran- 
dom, fitness dependent expected degree sequence [27, 28], 
which in this case is not known a priori and will be 
determined by the fitness distribution at the stationary 



^ log ^ = A^ 

ZT^ 



(9) 



If z remains finite as A^ — > oo, or in other words if zN — > 
cx), then the trivial solution is r — > 0. On the other 
hand, we find t 7^ if zN remains finite as A^ ^ 00. To 
obtain the value of t in this case, note that for a nonzero 
solution the term 1 /r in the above expression is finite and 
negligible for A^ large enough. Multiplying both sides by 
z yields 



^ log ^ = zN 



whose solution is r 



zN ' 



where 



(10) 



is the Pro- 



ductLog function defined as the solution of ^e"* = x. 
Putting these results together, we have 



(l){zN) 
zN 



1 



zN ^0 
y/(j){d)/d zN = d 



(11) 







zN 



As for random graphs, we find a marked transition 
as the scaling of z changes from N"^ to more rapidly 
decaying. This suggests an analogous underlying per- 
colation transition. As we show below, this is indeed 
the case. We can therefore still refer to the subcritical, 
sparse and dense regimes. Note that as N —>■ 00 we have 
/(x, y) = zxy in the sparse and subcritical regimes since 
zxy < z <^ I, which implies that we can neglect zxy 
in the denominator of eq.(l7|). Therefore the expression 
{f{x,m)) ~ f{x,{m)) is exact. On the other hand, in 
the dense regime we have r ^ 0, which again implies 
the same expression since q(m) becomes the Dirac delta 
function d{m). Therefore our trick to use the above ex- 
pression turns out to be exact in all regimes for N 00. 
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FIG. 1: Stationary fitness distribution and threshold. 

Main panel: cumulative density function p>{x) in log-linear 
axes. From right to left, z — 0.01, z = 0.1, z = 1, z = 10, 
z = 100, z = 1000 (iV = 5000). Inset: log-log plot of r{zN). 
Solid lines: theoretical curves, points: simulation results. 



Results and discussion 

The main panel of figd] shows the cumulative density 
function (CDF) of the fitness p> (x), while the inset shows 
a plot of t{zN). The theoretical results are in excellent 
agreement with numerical simulations. As predicted by 
eq.®, p{x) is the superposition of a uniform distribu- 
tion and a power-law with exponent —1. For z <^ 1 we 
have f(x,y) « zxy and p{x) oc x"^ for x > t. This 
purely power-law behaviour, that becomes exact in the 
sparse regime z = d/N for iV — > cx), results in a log- 
arithmic CDF looking like a straight line in log-linear 
axes. Note that, despite the value of the exponent, the 
presence of a nonzero lower threshold ensures that p{x) 
is normalizable. This mechanism may provide a natural 
explanation for the onset of Pareto distributions with a 
finite minimum value in real systems. By contrast, for 
large z the uniform part is nonvanishing and p{x) devi- 
ates from the purely power-law behaviour. The decay 
of p{x) for a; > T is a completely novel outcome of the 
extremal dynamics due to the feedback with the topol- 
ogy: now the fittest species at a given time is also the 
most likely to be connected to the least fit species and to 
mutate at the following timestep. Being more connected 
also means being more subject to changes. This enriches 
the coexistence patterns displayed on static networks. 

We now check the conjectured percolation transition. 
For different system sizes, we find that the cluster size 
distribution P{s) displays power-law tails when the con- 
trol parameter d = zN approaches a critical value dc — 
1.32 ±0.05 (corresponding to Zc = dc/N), suggesting the 
onset of a second-order percolation-like phase transition. 
As shown in figure [H P{s) cx with 7 = 2.45 ± 0.05 
at the phase transition. Fig. [3] shows that the aver- 
age fraction of nodes in the largest component remains 
negligible for d < d^, whereas it takes increasing finite 
values above d,,. As an additional check, following the 
method adopted in ref.fs^, we have plotted the average 





-hd = 1.32 




□ d = 4 


' -t 
■ '-' + 


Od = 0.1 










□ ^^t. 








CD '^m 




□ 





FIG. 2: Cluster size distribution. Far from the critical 
threshold {d = 0.1 and d = 4), P{s) is well peaked. At 
dc = 1.32, P{s) oc s"^ with 7 = 2.45 =b 0.05. Here iV = 3200. 
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FIG. 3: Behaviour at the percolation threshold. Main 
panel: the fraction of nodes in the giant component for dif- 
ferent network sizes as a function of d. Inset: the non-giant 
component average size as a function of d for A'^ = 6400. 



size fraction of non-giant components, which diverges (in 
the infinite volume limit) when P(s) decays algebraically 
as reported in the inset of fig. [3] 

Even if one of the most studied properties of the BS 
model on regular lattices is the statistics of avalanches 
characterizing the SOC behaviour [l^, we do not con- 
sider it here. This is because, as shown in ref . [3lj| . when 
considering long-range (22] instead of nearest-neigbour 
connections, it can lead to a wrong assessment of the 
SOC state, which is put into question by the absence 
of spatial correlations even in the case that avalanches 
are power-law distributed. Rather, we further charac- 
terize the topology at the stationary state by considering 
the degree distribution P{k) and the degree correlations. 
Using eq.([ni), we find that the average degree k{x) of a 
vertex with fitness x is 

, , , 2 , 1 + zx zx — ln(l -f zx) 
k{x) = —\n— + ^ (12) 

•^rj-^ I I ^ ^ 



1 + ZTX 



ZTX 



Similarly, through eq.(|4]) we can determine the analyt- 
ical expression for the cumulative degree distribution 
Py{k). As shown in figlH k{x) is linear for small z 
since /(x, z) ~ zxy, while for large z it saturates to 
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FIG. 4: Fitness— dependence and cumulative distribu- 
tion of the degrees. Left: k(x) {N — 5000; from right to 
left, z = 0.01, z = 0.1, z = 1, z ^ 10, z ^ 100, z = 1000). 
Right: P>{k) (same parameter values, inverse order from left 
to right). Solid hnes: theoretical curves, points: simulation 
results. 



the maximum value kmax ~ This implies that in 

the sparse regime P(fc) mimics p(x) and is characterized 
by the threshold value fc(r) and by a power-law decay 
P{k) oc fc^^ above it (see figH]). Note that here r re- 
mains finite even if P{k) oc k~'^ with 7 < 3, in strik- 
ing contrast with what obtained on static scale-free net- 
works [23, H^, HI] ■ Differently, for large z the saturation 
k — > kmax translates into a cut-off that makes P{k) de- 
viate from thepure power-law behavior for k > k{T). As 
shown in refs.piH^I, this saturation determines anticorre- 



lation between the degrees of neighboring vertices (disas- 
sortativity) and a hierarchy of degree-dependent cluster- 
ing coefficients as observed in many real-world networks 
(this is not shown here for brevity). As TV ^ cx3, these 
correlations vanish in the sparse regime (t > 0), while 
they survive in the dense regime (r — > 0). Structural 
correlations and a nonzero threshold r are then mutually 
excluding in this model, which is another interesting ef- 
fect of the feedback we have introduced. We finally note 
that with a different choice of /, one can have any ex- 
ponent for the power-law distribution of fitnesses, and 
therefore for the degree distribution as well. This makes 
our model completely flexible in order to reproduce any 
desired topological property. 

Our results represent a first step into the unexplored 
domain of systems with generic self-organized coupling 
between dynamics and topology. A huge class of such 
processes needs to be studied in the future, to further 
understand the unexpected effects of this coupling. 
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